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1. Introduction 

Systems containing more than two mesons are of phenomenological interest in a number of 
areas from heavy ion collisions at RHIC, to the equation of state of neutron stars. In the last few 
years, there has been a concerted effort by the NPLQCD collaboration to study such systems from 
first principles using lattice QCD [1, 2, 3]. Much progress has been made and these studies enabled 
the first calculation of a three hadron interaction, that between three like-charged pions. 

On the theoretical side, a number of developments have also occurred in the last few years. 
First, previous results for the volume dependence of the ground state energies of n boson systems 
[4, 5, 6] have been extended to the case of mixed species systems involving pions and kaons [8]. 
Secondly, a new algorithm has been presented [9] that recursively generates the quark propagator 
contractions for complicated systems from those of smaller systems. In principle, this method 
allows the calculation of systems with very large numbers of mesons (surmounting the limit of 
n = 12 particles in previous calculations) and of systems involving many different types of bosons. 
In this contribution, recent aspects of this progress are highlighted. 

2. Multi-Meson Interactions 

It has long been known how to exploit the volume dependence of the eigen-energies of two 
hadron systems to extract infinite volume scattering phase shifts [10] provided that the effective 
range of the interaction, r is small compared to the spatial extent of the lattice volume (since 
r ~ m„ for most interactions, this constraint is equivalent to m n L 3> 1). In recent works, this has 
been extended to systems involving n > 2 bosons [4, 5, 6] and n = 3 fermions [7] in the situation 
where the relevant scattering length, a, is small compared to the spatial extent of the lattice. By 
performing a perturbative calculation involving two- and three- body hadronic interactions in a 
finite volume, the ground state energy of a system of n bosons was computed in Refs. [4, 5, 6]. The 
shift in energy of n bosons of mass M from the non-interacting system is 





) 



4 



J A - bJ 2 J + (4 + n - n 2 ) J 2 +4(27 - I5n + n 2 )J X 




+ "C 3 



192 a 5 
M7i 3 L 7 



M 3 L 7 



6na 



(n + 3) J 




(2.1) 



where the parameter a is related to the scattering length, a, and the effective range, r, by 
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The geometric constants, J^, , J?T, ££ , ,% \, that enter into Eq. (2.1) are defined in Ref. [6] and 
"C m are the binomial coefficients. The three-body contribution to the energy-shift given in Eq. (2. 1) 
is represented by the parameter T] 3 (see Ref. [6]). 

Using lattice QCD calculations of these energy shifts, the parameters a and r\ 3 can be ex- 
tracted (this method was used in Refs. [1, 2, 3]). To determine the energy shifts, the multi-meson 
correlation functions (specifying to the multi-pion system) 

c »(0 « ((l>~ (*>')) (^ + (0,0)) ^ , (2.3) 
are calculated. On a lattice of infinite temporal extent, 1 the combination 

G - w = riff t ^ ; ^ n)e ~ AEnt > (2 - 4) 

allows for an extraction of the ground-state energy shift, AE n , which can then be used as input 
into Eq. (2.1) to extract the scattering and interaction parameters. To compute the (n\) 2 Wick 
contractions in Eq. (2.3), the correlation function can be written as 

c n (t) oc ( ( tjtit] T ) - e**-**-&-. e m . Ml . 4n _ n (n)* (n)g .. (n)£ , 

n = £ 5(x,r,0,0) S f (x,*;0,0) , (2.5) 

X 

where S(x,t;0,0) is a light-quark propagator. The object (block) H is a. N x N bosonic time- 
dependent matrix where N = 12 (/V = N$ x Nc with Ns = 4-spin and Nc = 3-color), and r\ a is a 
twelve component Grassmann variable. Further simplifications are possible resulting in the corre- 
lation functions being written in terms of traces of powers of n. As an example, the contractions 
for the 3-7T + system give 

C 3 (0 ~ tr [n] 3 - 3 tr [n 2 ] tr [n] +2 tr [n 3 ] , (2.6) 
where the traces are over color and spin indices. 



3. Contractions for Large (N > 12) Systems of Mesons 

To extend these types of calculations to systems of n > 12 mesons, or to study systems of many 
different types of mesons, different methods of performing the contractions of quark fields are 
required. While better than the naive factorial construction (which scales as n\ ), the construction 
used in Section 2 scales poorly to large numbers of mesons, behaving at best as n\ 1 / 2 (provided the 
matrix IT is generalized to give a nonzero result). In Ref. [9], a recursive method for performing 
these contractions was developed that allows the extension of the study of meson systems to larger 
n and also greatly simplifies the contractions required for systems of many different species of 
mesons. We will outline the construction by considering the recursive approach to the contractions 
for a single species of meson before reporting the general case. 

'Effects of temporal (anti-)periodicity are discussed in Ref. [3]. 
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By rescaling the correlation functions of the single-species, single-source system considered 
previously as 

C nJC+ (0 = (-)"«! (R n ) , (3.1) 

(the angle brackets denote a trace over spin and color indices) it is straightforward to show that the 
objects, R n that are defined implicitly in the above equation satisfy an ascending recursion 

R„+i = (R„)U - nR n U. (3.2) 

with the initial condition that R\ = IT as defined in Eq. (2.5). To see how this works we explicitly 
construct the first few terms: 

r 2 = { Ri )n - Ri n = ( n ) n - n 2 
<* 2 ) = <n) 2 - <n 2 > , 

r 3 = { r 2 ) n - 2 r 2 n = ( n ) 2 n - { n 2 ) n - 2 { n ) n 2 + 2 n 3 
( r 3 ) = ( n ) 3 - 3 ( n 2 ) ( n ) - 2 < n 3 > , (3.3) 

in agreement with Eq. (2.6). Descending recursions also exist. 

A correlation function for a system composed of nij mesons of the i th species from the / h 
source at (y/,0), where < i < k and < j < m, is of the form 

/ \ 111 / \n\ m / \ n k \ / \n km \ 

K( yi ,0)J ... K(y m ,0)J ... K(yi,0)J ... K t (y m ,0)J \ , (3.4) 

where JV{ = Y,j nij is the total number of mesons of species i, and the subscript in C n (t) labels the 
number of each species from each source, 
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(3.5) 



\n k \ n k2 ... n k „ 



The £/i(y,t) denotes a quark-level interpolating operator £/ m (x,t) = q m (x,t) 75 u(x,t), and it can 
be shown that [9] 

cAt) = (n^ ! ) I^UinPijnr^ = ( )■ i{h iil] ^"'^ {Tn ) ( 3.6) 

where the r\ are m x A^-component Grassmann variables, and the Pij are N xN dimensional matri- 
ces, where N = m x N, which are generalizations of the n defined in Eq. (2.5) with an additional 
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species index, i. They are denned as 
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(3.7) 



The ( A, ) ah axe N x N dimensional matrices, one for each flavor, i, and pair of source indices, a 
and b. JV = Ya 'M is tne tota l number of mesons in the system, with Jf < N. The T n defined 
implicitly in Eq. (3.6) satisfy the recursion relation 



T B +l rs = £ £ i r n+l„-ly ) Pij 
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(3.8) 



where 
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(3.9) 







and where the sole non-zero value is in the (/,j') th entry. Defining = ^ to be the number of 
mesons from the j th source, it is clear that the correlation function vanishes when %j > N for any 
source j. 

These recursion relations (and those for the simpler systems also studied in Ref. [9]) allow 
for the calculation of arbitrarily large systems of mesons. As formulated above, the systems are 
restricted to contain quarks of one flavor but anti-quarks of any number of flavors 2 or vice versa. 
Importantly the above algorithm scales linearly with the total number of mesons in the system. 

This recursive construction has been implemented and is currently being used to study a variety 
of complex meson systems that were previously not accessible. Generalisations to the case of 
baryons are straightforward but cumbersome. 



4. Mixed Species Meson Systems 

In the case of a mixed system comprising n mesons of one type and m mesons of a second type, 
the results summarized in Section 2 have been generalized in Ref. [8], working to ff(\/L ). The 
energy shift of the interacting system from the free system depends on three two-body interaction 
parameters and four three-body interaction parameters. As might be expected, the full form of the 
energy shift is lengthy and we refer the interested reader to the original paper. 

2 This restriction is not necessary, but relaxing it results in much more complex sets of recursions. 
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These systems are currently under numerical study and preliminary results using L 3 x T = 
20 3 x 128 anisotropic clover lattices at a pion mass of 390 MeV generated by the Hadron Spectrum 
Collaboration are discussed herein. 

In order to extract the energies of the multi meson systems, we study the behavior of two point 
correlation functions 

c N ,M(t) = ( ( E ir (x, t)) N ( £k (x, t)) M (k + (o, t)) N (k+ (o, t)) M ) , (4. 1) 

where the pion and kaon interpolating operators are defined in terms of quark fields as n + (x,t) = 
u(x,t)ysd(x,t) and K + (x,t) = u(x,t)ysd(x,t), respectively. 

At large Euclidean times, <C t <C T , these correlation functions decay exponentially with the 
ground state energies, Emm, of the system. These energies could then be used to constrain the two- 
and three- body interactions using the results of Ref. [8]. However, satisfying the constraint <C 
t <T, requires taking the limit of vanishing temperature and is not typically practical in current 
lattice calculations. The presence of thermal fluctuations that vanish exponentially as T — > °°, leads, 
at finite temporal extent, to significantly more complex behaviour of the correlation functions in 
Eq. (4.1). The general form for these correlation functions can be straightforwardly deduced from 
Hamiltonian evolution at finite temperature and is given by the following 

i M N 

CNM(t) = ? £ £ Z^M-m^^+E^T/l CQsh U EN _ nM _ m -E njm ){t-T/2)) 

N M / \ . 

+-Zfl e ( E ^r-)V2 §N2l s Mak + . . . (4 . 2) 

Z 2 i 2 

where the last term is only present when N and M are both even numbers. The ellipsis denotes 
excited state contributions. The subleading terms shown in the summation correspond to the various 
ways one can connect the A^-pion, M-kaon source and sink, allowing for propagation of subsets of 
hadrons around the temporal boundary (various of the are simply related). 

With the large number of parameters that occur in Eq. (4.2), fitting the correlation functions be- 
comes a complicated endeavour. Our current analysis strategy employs variable projection (VarPro 
[12, 13]) to eliminate the linear parameters, the Z^ m "' M ~ m , from the minimisation procedure. As ,/V 
and M increase, we also employ the method of Bayesian priors [14] by providing suitable physical 
starting points for the value of the energies in the correlation function. A summary plot of the 
energy extractions is shown in Fig. 1 . 

Analysis of these correlation functions is ongoing and preliminary results are very encour- 
aging. The data provide statistically significant extractions of most of the two- and three- body 
interactions that enter the mixed-species systems. It additionally allows us to study the relation- 
ships between the isospin and hypercharge densities and the corresponding chemical potentials. 

5. Summary 

Systems of large numbers of mesons provide a novel form of QCD matter that is amenable 
to study using lattice QCD. A rich phase structure is expected and the tools and methodologies 
presented herein will allow a close study of these systems. These techniques also represent an 
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Figure 1: Energies of systems of N pions and M kaons extracted from our calculations. The vertical extent 
of the shaded boxes corresponds the uncertainty in the determination of the energies. 

important step towards the complexity frontier that must be confronted if lattice QCD is to be 
applied to the physics of nuclei. 
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